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Abstract.1  Proper  control  of  the  numerical-dissipation/filter  to  accurately  resolve  all  relevant 
multiscales  of  complex  flow  problems  while  still  maintaining  nonlinear  stability  and  efficiency  for 
long-time  numerical  integrations  poses  a great  challenge  to  the  design  of  numerical  methods.  The 
required  type  and  amount  of  numerical-dissipation/filter  are  not  only  physical  problem  depen- 
dent. but  also  vary  from  one  flow  region  to  another.  An  approach  for  the  automatic  detection  of 
different  flow  features  as  distinct  sensors  to  signal  the  appropriate  type  and  amount  of  numerical- 
dissipation/filter  for  non-dissipative  high  order  schemes  is  proposed.  These  scheme-independent 
sensors  are  capable  of  distinguishing  shocks/shears,  turbulent  fluctuations  and  spurious  high- 
frequency  oscillations.  In  addition,  these  sensors  are  readily  available  as  an  improvement  over 
existing  grid  adaptation  indicators.  The  same  shock/shear  detector  that  is  designed  to  switch  on 
the  shock/shear  numerical  dissipation  can  be  used  to  switch  off  the  entropy  splitting  form  of  the 
inviscid  flux  derivative  [31]  in  the  vicinity  the  discontinuous  regions  to  further  improve  nonlinear 
stability  and  minimize  the  use  of  numerical  dissipation.  The  rest  of  the  sensors  in  conjunction 
with  the  local  flow  speed  and  Reynolds  number  can  also  be  used  to  adaptively  determine  the 
appropriate  entropy  splitting  parameter  for  each  flow  type/region.  The  goal  of  this  paper  is  to 
further  improve  nonlinear  stability,  accuracy  and  efficiency  of  long-time  numerical  integration  of 
complex  shock/turbulence/acoustics  interactions  and  numerical  combustion.  The  minimization 
of  employing  very  fine  grids  to  overcome  the  production  of  spurious  numerical  solution  and/or 
instability  due  to  under-resolved  grids  is  also  sought  [29,  6]. 


1.  Motivation  and  Objective 

When  employing  finite  time  steps  and  finite  grid  spacings  in  the  long-time  integration  of  mul- 
tiscale complex  nonlinear  fluid  flows,  nonlinear  instability  commonly  occurs  and  the  use  of  a 
numerical-dissipation/filter  is  unavoidable.  Aside  from  acting  as  a post-processor  step,  most 
filters  serve  as  some  form  of  numerical  dissipation.  Without  loss  of  generality,  “numerical- 
dissipation/filter”  is,  hereafter,  referred  to  as  “numerical  dissipation”  unless  otherwise  stated. 
The  required  type  and  amount  of  numerical  dissipation  can  vary  drastically  from  one  flow  region 
to  another  (within  the  same  problem)  and  are  highly  problem  dependent.  This  is  particularly 
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true  for  unsteady  high-speed  turbulence  and/or  combustion  problems.  For  this  type  of  problem, 
it  is  of  paramount  importance  to  have  proper  control  on  the  type  and  amount  of  numerical 
dissipation  in  regions  where  it  is  needed  but  nowhere  else.  Inappropriate  type  and/or  amount 
can  be  detrimental  to  the  integrity  of  the  computed  solution,  especially  for  problems  with  no 
known  analytical  solution  and/or  experimental  data  for  comparison. 

The  linear  and  nonlinear  numerical  dissipation  presently  available  is  either  built  into  the 
numerical  scheme  or  added  on  to  the  existing  scheme.  The  built-in  numerical  dissipation  schemes 
are,  e.g.,  upwind,  flux  corrected  transport  (FCT),  total  variation  diminishing  (TVD),  essentially 
non-oscillatory  (ENO),  weighted  ENO  (WENO),  and  hybrid  schemes  (e.g.,  those  that  switch 
between  spectral  and  high-order  shock-capturing  schemes).  The  built-in  nonlinear  numerical 
dissipation  in  TVD,  ENO  and  WENO  schemes  was  designed  to  capture  accurately  discontinuities 
and  high  gradient  flows  while  hoping  to  maintain  the  order  of  accuracy  of  the  scheme  away  from 
discontinuities. 

There  exist  different  specialized  linear  and  nonlinear  filters  to  post  process  the  numerical 
solution  after  the  completion  of  a full  time  step  of  the  numerical  integration.  Since  they  arc 
post  processors,  the  physical  viscosity,  if  it  exists,  is  taken  into  consideration.  The  main  design 
principle  of  linear  filters  is  to  improve  nonlinear  stability,  under-resolved  grids  [6]  and  de-aliasing 
of  smooth  flows,  while  the  design  principle  of  nonlinear  filters  is  to  improve  the  nonlinear  stability 
and  accuracy  near  discontinuities  as  well.  See,  for  example,  the  present  proceedings  and  [9,  6, 
7,  27]  for  forms  of  linear  filters,  and  see  [30,  31,  21]  for  forms  of  nonlinear  filters.  For  direct 
numerical  simulation  (DNS)  and  large  eddy  simulation  (LES),  there  are  additional  variants  of 
the  filter  approach. 

For  the  last  decade,  CPU  intensive  high  order  schemes  with  built-in  nonlinear-  dissipation 
have  been  gaining  in  popularity  in  DNS  and  LES  for  long-time  integration  of  shock-turbulence 
interactions.  Unfortunately,  these  built-in  dissipations  cannot  clearly  distinguish  shocks/shears, 
from  turbulent  fluctuations  and/or  spurious  high-frequency  oscillations.  In  [30,  31,  21),  it  was 
shown  that  these  built-in  numerical  dissipation  are  more  dissipative  and  less  accurate  than  the 
nonlinear  filter  approach  of  [30,  31,  21]  with  a similar  order  of  accuracy.  It  was  also  shown 
that  these  nonlinear  filters  can  also  suppress  spurious  high-frequency  oscillations.  However,  a 
subsequent  study  of  Sjogreen  & Yee  [23]  showed  that  the  high  order  linear  filter  can  sustain 
longer  time  integration  more  accurately  than  the  nonlinear  filter  for  low  speed  smooth  flows.  In 
other  words,  for  the  numerical  examples  that  were  examined  in  [23],  the  high  order  linear  filter 
can  remove  high-frequency  oscillations  producing  nonlinear  instability  better  than  the  second- 
order  nonlinear  filter.  Higher  than  third-order  nonlinear  filters  might  be  able  to  improve  their 
performance  or  might  outperform  the  high  order  linear  filters  at  the  expense  of  more  CPU  time 
and  added  complexity  near  the  computational  boundaries.  These  findings  prompted  the  design 
of  switching  on  and  off  or  blending  of  different  filters  to  obtain  the  optimal  accuracy  of  high  order 
spatial  difference  operator  as  proposed  in  Yee  et  al.  and  Sjogreen  & Yee  [31,  21).  The  missing  link 
of  what  was  proposed  in  [31,  21]  is  an  efficient,  automated  and  reliable  set  of  appropriate  sensors 
that  are  capable  of  distinguishing  shocks/shears,  from  turbulent  fluctuations  and/or  spurious 
high-frequency  oscillations  for  a full  spectrum  of  flow  speeds  and  Reynolds  numbers. 

The  present  paper  is  a sequel  to  [30,  31,  21).  The  objective  here  is  to  propose  an  adaptive 
procedure  employing  appropriate  sensors  to  switch  on  the  desired  numerical  dissipation  where 
needed  and  leave  the  rest  of  the  region  free  of  numerical  dissipation  contamination  while  at 
the  same  time  improving  nonlinear  stability  of  the  entire  numerical  process.  In  addition,  the 
minimization  of  employing  very  fine  grids  to  overcome  the  production  of  spurious  numerical 
solutions  and/or  instability  due  to  under-resolved  grids  is  sought  [6].  It  was  shown  in  [17,  8,  31,  21] 
that  conditioning  the  governing  equations  via  entropy  splitting  of  the  inviscid  flux  derivatives  [31] 
can  improve  the  over  all  stability  of  the  numerical  approach  for  smooth  flows.  Therefore,  the  same 
shock/shear  detector  that  is  designed  to  switch  on  the  shock /shear  numerical  dissipation  can 
be  used  to  switch  off  the  entropy  splitting  form  of  the  inviscid  flux  derivative  in  the  vicinity  the 
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discontinuous  regions  to  further  improve  nonlinear  stability  and  minimize  the  use  of  numerical 
dissipation.  The  rest  of  the  sensors,  in  conjunction  with  the  local  flow  speed  and  Reynolds 
number,  can  also  be  used  to  adaptively  determine  the  appropriate  entropy  splitting  parameter 
for  each  flow  type/region.  These  sensors  are  readily  available  as  an  improvement  over  existing 
grid  adaptation  indicators.  If  applied  correctly,  the  proposed  adaptive  numerical  dissipation 
control  is  scheme  independent,  and  can  be  a stand  alone  option  for  many  of  the  favorite  schemes 
used  in  the  literature.  Although  the  proposed  adaptive  procedure  is  scheme  independent,  we 
prefer  a complete  treatment  of  the  numerical  approach  in  the  following  framework: 

1.  For  stability  considerations,  condition  the  governing  equations  before  the  application  of  the 
appropriate  numerical  scheme. 

2.  For  consistency,  compatible  schemes  that  posses  stability  properties  similar  to  those  of  the 
discrete  analogue  of  the  continuum  are  preferred. 

3.  For  the  minimization  of  numerical  dissipation  contamination,  efficient  and  adaptive  nu- 
merical dissipation  control  to  further  improve  nonlinear  stability  and  accuracy  should  be 
used. 

4.  For  practical  considerations,  the  numerical  approach  should  be  efficient  and  applicable  to 
general  geometries.  An  efficient  and  reliable  dynamic  grid  adaptation  should  be  used  if 
necessary.  Note  that  the  computation  of  the  illustrative  examples  used  such  a numerical 
approach  [30,  31,  28,  21,  8,  16]. 

A brief  summary  of  (l)-(3)  is  discussed  in  the  next  three  sections.  Some  representative  examples 
to  illustrate  the  performance  of  the  approach  are  given  in  Section  5. 

2.  Conditioning  of  the  Governing  Equations 

Traditionally,  conditioning  the  governing  partial  differential  equations  (PDEs)  usually  referred 
to  rewriting  the  governing  equations  in  an  equivalent  set  of  PDEs  in  order  to  prove  the  stability 
and/or  wcll-posedness  of  the  PDEs.  When  numerical  methods  are  used  to  solve  PDEs  that 
are  nonlinear,  it  is  well-known  that  different  equivalent  forms  of  the  governing  equations  might, 
exhibit  different  numerical  stability,  accuracy  and/or  spurious  computed  solutions,  even  for 
problems  containing  no  shock/shear  discontinuities.  There  are  many  conditioned  forms  of  the 
governing  equations  proposed  in  the  literature.  Different  conditioned  forms  of  the  convection 
fluxes  and  the  viscous  fluxes  have  been  proposed  for  the  incompressible  and  compressible  Navier- 
Stokes  equations.  Here  we  mention  a few  which  are  precursors  of  the  so  called  entropy  splitting 
of  the  compressible  Euler  equations  [31].  If  a mcthod-of-lincs  approach  is  used  to  discretize  these 
equations,  the  entropy  splitting  reduces  to  the  splitting  of  the  convection  flux  derivatives. 

The  splitting  of  the  convection  terms  (for  both  the  compressible  and  incompressible  Navier- 
Stokes  equations)  into  a conservative  part  and  a non-conservative  part  has  been  known  for  a 
long  time.  In  the  DNS,  LES  and  atmospheric  science  simulation  literature,  it  is  referred  to  as  the 
skew-symmetric  form  of  the  momentum  equations  [1,  14,  2,  35].  It  consists  of  the  mean  average  of 
the  conservative  and  non-conservative  (convective  form  [35])  part  of  the  momentum  equations. 
The  spatial  difference  operator  is  then  applied  to  the  split  form.  From  the  numerical  analysis 
standpoint,  the  Hirt  and  Zalesak’s  ZIP  scheme  [12,  34]  is  equivalent  to  applying  central  schemes 
to  the  non-conservative  momentum  equations  (convective  form  of  the  momentum  equations). 
MacCormack  [15]  proposed  the  use  of  the  skew-symmetric  form  for  problems  other  than  DNS 
and  LES.  Harten  [10]  and  Tadmor  [26]  discussed  the  symmetric  form  of  the  Euler  equations 
and  skew-adjoint  form  of  hyperbolic  conservation  laws,  respectively.  Although  the  derivation  in 
these  works  is  different,  the  ultimate  goal  of  using  the  split  form  is  almost  identical.  This  goal  is 
to  improve  nonlinear  stability,  minimize  spurious  high-frequency  oscillations  and  robustness  of 
the  numerical  computations.  See  [33,  23]  for  a historical  consolidation  of  these  approaches.  The 
canonical  splitting  used  by  Olsson  & Oliger  [17]  is  a mathematical  tool  to  prove  the  existence 
of  a generalized  energy  estimate  for  a symmetrizable  system  of  conservation  laws.  For  the  ther- 
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mally  perfect  gas  compressible  Euler  equations,  the  transformation  consists  of  a convex  entropy 
function  that  satisfies  a mathematical  entropy  condition.  The  mathematical  entropy  function, 
in  this  case,  can  be  a function  of  the  physical  entropy.  Therefore,  the  resulting  splitting  was 
referred  to  as  entropy  splitting  for  ease  of  reference  by  Yee  et  al.  [31].  The  entropy  splitting 
can  be  viewed  as  the  more  general  form  which  provides  L2  stability  proof  of  the  nonlinear  Euler 
equations  with  physical  boundary  conditions  included. 

Consider  a general  nonlinear  system  of  conservation  laws, 

Ut  + F(U)x  = 0,  a < x <b.  (2.1) 

If  the  conservation  law  has  an  entropy  function,  it  can  be  transformed  into  a symmetric  con- 
servation law  in  terms  of  a new  variable  W.  The  change  of  variables  dU /dW  is  symmetric  and 
positive  definite,  and  the  new  Jacobian  dF/dW  is  symmetric.  If  furthermore  U and  F are  ho- 
mogeneous in  W of  degree  ft,  which  is  the  case  for  the  thermally  perfect  gas  Euler  equations  for 
any  ft  5^  — 1,  the  formulas  become  simple.  In  that  case  we  insert  the  change  of  variables  into  the 
conservation  law,  and  define  the  split  form  of  the  flux  derivative  [17] 

u>  + rhF*+TbFwWr=0'  (2-2) 

with  ft  a splitting  parameter  (ft  = 00  recovers  the  original  conservative  form).  Here  ft  — 1 and, 

for  a perfect  gas,  ft  > 0 or  ft  < The  theory  only  gives  the  range  of  ft  and  does  not  give  any 
guidelines  on  how  to  choose  ft  for  the  particular  flow.  The  vectors  Fw  and  W can  be  cast  as 
functions  of  the  primitive  variables  and  ft.  FVom  the  study  of  [31],  ft  is  highly  problem  dependent. 
Multiplying  the  above  equation  by  W (where  (p,q)  denotes  the  scalar  or  inner  product  of  the 
vectors  p and  q) , gives 

-(1  + ft)(W,  Ut ) = ft(W,  Fx)  + (JY,  FWWX)  = ft(W,  Fx)  + (FwW,  Wx). 

Integration  by  parts  in  space  gives 

(l  + ft)(W,Ut)  = -[WTFwW]ba 
from  which  we  obtain  the  estimate 

£(W,  UwW)  = (WL,  UwW)  + (W,  (UwW)t)  = (Ut,  W)  + ft(W,  Ut ) = 

(1  +ft){W,Vt)  = -[W'rFwW]ba. 

In  order  to  have  an  energy  estimate,  the  boundary  term  \WT FwW}^  should  be  of  the  sign 
that  makes  the  time  derivative  of  the  norm  negative.  For  stability  the  entropy  norm  {W,  IJ^W) 
should  be  bounded. 

3.  Discrete  Analogue  of  the  Continuum 

For  ease  of  reference,  “scheme”  or  more  precisely  “interior  scheme”  here  generally  refers  to  spatial 
difference  schemes  for  the  interior  grid  points  of  the  computational  domain,  whereas  “boundary 
scheme”  is  the  numerical  boundary  difference  operators  for  grid  points  near  the  boundaries. 
However,  without  loss  of  generality,  we  also  adopt  the  conventional  terminology  of  denoting 
“scheme”  interchangeably  as  either  the  “combined  interior  and  boundary  scheme”  or  just  the 
“interior  scheme”  within  the  context  of  the  discussion.  Before  1994,  rigorous  stability  estimates 
for  accurate  and  appropriate  boundary  schemes  associated  with  fourth-order  or  higher  spatial 
interior  schemes  were  the  major  stumbling  block  in  the  theoretical  development  of  combined 
interior  and  boundary  schemes  for  nonlinear  systems  of  conservation  laws.  Olsson  [18]  proved 
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that  an  energy  estimate  can  be  established  for  second-order  central  schemes.  To  obtain  a rigorous 
energy  estimate  for  high  order  central  schemes,  one  must  apply  the  scheme  to  the  split  form 
of  the  inviscid  governing  equation.  A discrete  analogue  of  the  continuum  using  a semi  discrete 
approach  can  be  written  as 

^‘-ihDFiB^rhFMUt)DWl‘  (31) 

Here,  D is  a difference  operator,  having  the  summation  by  parts  (SBP)  property  [18,  25].  The 
estimate 

i(W,  UwW)h  = —WJfw{Wj)Wj  + W[Fw(W1)Wl  (3.2) 

at 

in  the  discrete  scalar  product  follows  in  the  same  way  as  for  the  PDE  with  indices  1 and  J 
the  end  points  of  the  computational  domain,  and  h the  grid  spacing.  Here  the  SBP  satisfying 
difference  operator,  for  example,  consists  of  central  difference  interior  operators  of  even  order 
together  with  the  corresponding  numerical  boundary  operators  that  obey  the  discrete  energy 
estimate.  See  Olsson  and  Strand  for  forms  of  the  SBP  boundary  operators  [18,  25]. 

For  the  full  discretization  of  the  problem,  we  should  discretize  in  time  in  such  a way  that 
the  discrete  energy  estimate  also  holds.  The  obvious  solution  would  be  to  discretize  in  time  in  a 
skew  symmetric  way,  in  a manner  similar  to  the  spatial  discretization,  e.g., 

J~oD<U?  + ^fw{Wf)DlW]  = -^DFiU")  - - ±-pFw(W?)DW?.  (3.3) 

However,  it  turns  out  that  the  SBP  property  of  the  time  difference  quotient  leads  to  a problem 
which  is  coupled  implicitly  in  the  time  direction.  To  solve  it  we  have  to  solve  a nonlinear  system 
of  equations  for  all  time  levels  in  the  same  system,  leading  to  an  impractically  large  computa- 
tional effort.  Furthermore,  numerical  experiments  shown  in  Sjogreen  & Yee  [23]  indicated  that  a 
bounded  L2  entropy  norm  ( W . U\yW )/,  does  not  necessarily  guarantee  a well  behaved  numerical 
solution  for  long-time  integrations.  In  other  words,  L 2 stability  does  not  necessarily  guarantee 
an  accurate  solution.  In  practical  computations,  the  classical  Runge-Kutta  time  discretizations 
using  the  method-of-lines  approach  (which  we  used  for  our  numerical  experiments)  works  well, 
but  we  have  not  been  able  to  prove  a time  discrete  entropy  estimate  for  this  method.  In  addition, 
numerical  experiments  shown  in  [23]  indicate  that  the  time  discrete  problem  does  not  have  a 
decreasing  entropy  norm  for  all  values  of  0.  Numerical  experiments  in  Yee  et  al.  [31,  19]  also 
indicate  the  wide  variations  of  the  0 value  for  a full  spectrum  of  flow  problems.  For  example, 
if  a constant  0 is  used  for  problems  containing  shock  waves,  a very  large  value  of  0 is  needed. 
Otherwise  diverge  solution  or  wrong  shock  location  and/or  shock  strength  are  obtained.  With 
t hese  findings,  employing  a constant  0 (within  the  allowable  range  of  0)  throughout  the  entire 
computational  domain  is  not  advisable  unless  the  flow  problem  is  a simple  smooth  flow.  Studies 
in  [31,  21]  indicate  that  the  split  form  of  the  inviscid  flux  derivatives  does  help  in  minimizing 
the  use  of  numerical  dissipation.  What  is  needed  is  adaptive  control  of  the  0 parameter  from 
one  flow  region  to  another  as  well  as  from  one  physical  problem  to  another. 

We  would  like  to  point  out  that  for  non-homogeneous  physical  boundary  conditions,  it  is 
important  to  impose  these  boundary  conditions  to  maintain  the  SBP  property.  For  a nonlin- 
ear symmetrizable  system  of  conservation  laws  that  contain  time-dependent  physical  boundary 
conditions,  an  extra  complication  arises,  especially  if  a multistage  Runge-Kutta  method  is  used. 
See  Johansson  [13]  for  a discussion. 

4.  Adaptive  Numerical  Dissipation  Control 

For  smooth  flows,  entropy  splitting  does  significantly  improve  the  stability  of  the  computa- 
tion, but  does  not  completely  remove  all  nonlinear  instabilities,  especially  for  long-time  wave 
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propagations,  computation  of  turbulent  statistics,  and  spurious  numerical  solutions  and  spurious 
oscillations  due  to  under-resolved  grids  [29,  6],  Furthermore,  the  entropy  splitting  was  developed 
for  smooth  solutions.  Additional  difficulties  from  shock  wave  formation  are  the  oscillatory  behav- 
ior of  the  centered  difference  interior  scheme,  and  the  fact  that  the  nonconservative  terms  might 
lead  to  inconsistent  behavior  at  shocks/shears  [31,  21,  23].  We  propose  to  enhance  the  above 
conditioning  of  the  equations  with  an  advanced  numerical  dissipation  model,  which  includes 
nonlinear  sensors  to  detect  shocks/shears  and  other  small  scale  features,  and  spurious  oscilla- 
tion instability  due  to  under-resolved  grids.  Furthermore,  we  will  use  the  detector  to  switch  off 
the  entropy  splitting  at  shocks/shears  and  adjust  the  entropy  splitting  parameter  with  the  aid  of 
local  Macli  number  and  Reynolds  number  in  smooth  region  as  discussed  earlier.  The  advanced 
numerical  dissipation  model  can  be  used:  (Option  I)  as  part  of  the  scheme,  (Option  II)  as  an 
adaptive  filter  control  after  the  completion  of  a full  time  step  of  the  numerical  integration  or 
(Option  III)  as  a combination  of  Options  I and  II  (e.g.,  high  order  nonlinear  dissipation  (with 
sensor  control)  using  Option  I and  nonlinear  filter  (with  a different  sensor  control)  using  Option 
II). 

Due  to  space  constraints,  we  concentrate  here  on  an  adaptive  procedure  that  can  distin- 
guish three  major  computed  flow  features  to  signal  the  correct  type  and  amount  of  numerical 
dissipation  needed  in  addition  to  controlling  the  entropy  splitting  parameter.  The  major  flow 
features  and  numerical  instability  are  (a)  shocks/shears  , (b)  turbulent  fluctuations  , and  (c) 
spurious  high-frequency  oscillations.  The  procedure  can  be  extended  if  additional  refinement  or 
classification  of  flow  types  and  the  required  type  of  numerical  dissipation  is  needed.  There  exist 
different  detection  mechanisms  in  the  literature  for  the  above  three  features.  These  detectors 
are  not  mutually  exclusive  and/or  are  too  expensive  for  practical  applications.  We  believe  that 
the  multiresolution  wavelet  approach  proposed  in  Sjogreen  & Yee  [21]  is  capable  of  detecting  all 
of  these  flow  features,  resulting  in  three  distinct  sensors.  If  chosen  properly,  one  multiresolution 
wavelet  basis  function  might  be  able  to  detect  all  three  features.  For  an  optimum  choice,  one 
might  have  to  use  more  than  one  type  of  wavelet  basis  functions  but  at  the  expense  of  an  in- 
crease in  CPU  requirements.  Some  incremental  studies  into  the  use  of  entropy  splitting  and  the 
application  of  these  sensors  were  illustrated  in  [31,  21,  22,  24,  23,  16].  Here,  additional  tools  are 
illustrated  together  with  highlights  of  some  examples.  An  adaptive  numerical  dissipation  model 
illustrating  some  of  the  proposed  ideas  will  be  described  in  the  next  section.  Full  implementation 
of  the  complete  treatment  of  the  numerical  approach  (including  grid  adaptation)  discussed  in 
Section  1 to  complex  realistic  multiscale  problems  is  in  progress. 


5.  Numerical  Examples 

This  section  illustrates  the  power  of  entropy  splitting,  the  difference  in  performance  of  linear 
and  nonlinear  (with  sensor  controls)  filters  and  the  combination  of  both  types  of  filters  with 
adaptive  sensor  controls.  We  use  the  same  notation  as  in  [30,  31,  22],  The  artificial  compression 
method  (ACM)  and  wavelet  filter  schemes  using  a second-order  nonlinear  filter  with  sixth-order 
spatial  central  interior  scheme  for  both  the  iuviscid  and  viscous  flux  derivatives  are  denoted  by 
ACM66  and  WAV66.  See  [30,  31,  22]  for  the  forms  of  these  filter  schemes.  The  same  scheme 
without  filters  is  denoted  by  CEN66.  The  scheme  using  the  fifth-order  WENO  for  the  inviscid 
flux  derivatives  and  sixth-order  central  for  viscous  flux  derivatives  is  denoted  by  WEN05.  Com- 
putations using  the  standard  fourth-order  Runge-Kutta  temporal  discretization  arc  indicated 
by  appending  the  letters  “RK4”  as  in  CEN66-RK4.  ACM66  and  WAV66  use  the  Roe’s  average 
state  and  the  van  Leer  limiter  for  the  nonlinear  numerical  dissipation  portion  of  the  filter.  The 
wavelet  decomposition  is  applied  in  density  and  pressure,  and  the  maximum  wavelet  coefficient 
of  the  two  components  is  used.  The  nonlinear  numerical  dissipation  is  switched  on  wherever  the 
wavelet  analysis  gives  a Lipschitz  exponent  [21]  less  than  0.5.  Increasing  this  number  will  reduce 
oscillations,  at  the  price  of  reduced  accuracy  (see  [21]  for  other  possibilities).  Computations 
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Isentropic  Vortex  Evolution 

(Horizontally  Convectlng  Vortex,  vortex  strength  @=5) 
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Periodic  BC  in  * & y 


Initial  Vortex,  Penalty  Contour* 


center  of  vortex 


Figure.  5.1.  Vortex  convection  problem  description. 


using  entropy  splitting  are  indicated  by  appending  the  letters  “ENT”  as  in  WAV66-RK4-ENT. 
Computations  using  an  eighth-order  linear  dissipation  filter  are  indicated  by  appending  the  let- 
ters “D8”  as  in  WAVGG-RK4-D8.  In  order  not  to  introduce  additional  notation,  inviscid  flow 
simulations  are  designated  by  the  same  notation,  with  the  viscous  terms  not  activated. 

5.1.  A 2-D  VORTEX  CONVECTION  MODEL  [30,  31,  21,  23] 

The  onset  of  nonlinear  instability  of  long-time  numerical  integration,  the  benefit  of  the  entropy 
splitting  and  the  difference  in  performance  of  linear  and  nonlinear  numerical  dissipations  in 
improving  nonlinear  stability  for  a horizontally  convecting  vortex  (see  Fig.  5.1)  can  be  found  in 
in  [30,  31,  21,  23].  We  summarize  the  result  here. 

To  show  the  onset  of  nonlinear  instability,  the  2-D  perfect  gas  compressible  Euler  equations 
are  approximated  by  CEN66-RK4  with  periodic  boundary  conditions  imposed  using  a 80  x 79 
grid  with  the  time  step  At  = 0.01.  Since  this  is  a pure  convection  problem,  the  vortex  should 
convect  without  any  distortion  if  the  numerical  scheme  is  highly  accurate  and  non-dissipative. 
Although  CEN66-RK4  is  linearly  stable,  the  test  problem  is  nonlinear  and  instability  eventually 
sets  in.  Almost  perfect  vortex  preservation  is  observed  for  up  to  5 periods  of  integrations  (5 
times  after  the  vortex  has  converted  back  to  the  same  position  - time  = 50).  Beyond  5 periods 
the  solution  becomes  oscillatory,  and  blows  up  before  the  completion  of  6 periods.  The  blow 
up  is  associated  with  an  increase  in  entropy  [23].  If  we  instead  use  the  entropy-split  form  of 
the  approximation  (CEN66-RK4-ENT)  with  a split  parameter  ft  = 1,  almost  perfect  vortex 
preservation  for  up  to  40  periods  can  be  obtained.  The  computation  remains  stable  for  up  to  67 
periods  before  it  breaks  down.  The  time  history  of  the  L2  entropy  norm  and  density  contours 
of  the  solution  after  5,  10,  30  and  G7  periods  using  CENGG-RK4-ENT  is  shown  in  Figs.  5.2  and 
5.3.  The  norm  is  decreasing,  although  the  instabilities  break  down  the  solution  after  67  periods. 
Using  the  second-order  nonlinear  filter  without  splitting  (ACM66-RK4  or  WAV66-RK4),  the 
solution  remains  stable  beyond  67  periods.  However,  the  numerical  solution  gradually  starts 
to  diffuse  after  20  periods.  If  we  use  the  nonlinear  filter  in  conjunction  with  entropy  splitting 
(ACM66-RK4-ENT  or  WAV66-RK4-ENT),  almost  perfect,  vortex  preservation  can  be  obtain  for 
up  to  120  periods  using  a split  parameter  ft  [31]. 
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Entropy  Norm  vs  Time,  CEN6S-RK4-ENT  |M 


Figure  5.2.  Entropy  norm  history  of  CEN66-ENT:  entropy  split  parameter  /3  = 1 and  80  x 79  grid. 


Figure  5.3.  Density  Contours  of  CEN66-ENT:  entropy  split  parameter  >9  = 1 and  80  x 79  grid. 


The  density  contours  of  the  solution  after  5,  10,  200  and  300  periods  for  the  unsplit  (/?  = oo) 
computation  using  the  eighth-order  linear  dissipation  (CEN66-RK4-D8)  are  shown  in  Fig.  5.4. 
The  linear  dissipation  term  (—dh7(D+D-)4Uj)  with  grid  spacing  h was  added  to  the  sixth-order 
base  scheme  to  discretize  the  convection  terms.  The  parameter  d is  a given  constant  ( d = 0.0002) 
and  is  scaled  with  the  spectral  radius  of  the  Jacobian  of  the  flux  function,  and  D+  and  are  the 
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Figure  5-4-  Density  contours  of  CEN66-RK4-D8  using  80  x 79  grid. 


forward  and  backward  difference  operators,  respectively.  This  numerical  dissipation  is  applied 
as  part  of  the  scheme  and  not  as  a post  processing  filter.  The  computation  can  be  run  for  300 
periods  without  breakdown.  However,  serious  degradation  of  accuracy  occurs  after  250  periods. 
For  this  particular  problem,  the  CEN66-RK4-D8  out  performed  the  ACM66-RK4-ENT  and 
WAV66-RK4-ENT  using  0=1.  Perhaps  using  a higher  than  third-order  nonlinear  filter  might 
improve  the  performance  of  the  ACM66-RK4-ENT  and  WAV66-RK4-ENT  at  the  expense  of  ari 
increase  in  CPU. 

5.2.  DNS  OF  3-D  COMPRESSIBLE  TURBULENT  CHANNEL  FLOW  [19] 

To  obtain  an  accurate  turbulent  statistics,  very  long-time  integration  and  highly  accurate  meth- 
ods are  required  for  this  DNS  computation.  This  numerical  example  illustrates  the  power  of 
entropy  splitting.  The  computation  employed  the  SBP-satisfying  boundary  difference  opera- 
tor with  the  fourth-order  central  interior  scheme  applied  to  the  split  form  of  the  inviscid  flux 
derivatives  CEN44-RK4-ENT  with  0 = 4.  The  fluid  mechanics  of  this  3-D  wall  bounded  isother- 
mal compressible  turbulent  channel  flow  has  been  studied  in  some  detail  by  Coleman  et,  al. 
[4],  They  showed  that  the  only  compressibility  effect  at  moderate  Mach  numbers  comes  from 
the  variation  of  fluid  properties  with  temperature.  They  used  a uniform  body  force  term  to 
drive  the  flow,  but  recommended  the  constant  pressure  gradient  approach  which  was  adapted 
by  Sandham  & Yee  [19].  For  simplicity,  the  fixed  pressure  gradient  problem  rather  than  the 
constant  mass  flow  problem  was  solved  using  fixed  fluid  properties.  Thus  only  the  wall  shear 
stress  and  mass  flow  rate  vary  during  the  simulation.  A Mach  number  of  0.1  is  chosen,  based 
on  friction  velocity  and  sound  speed  corresponding  to  the  fixed  wall  temperature.  Channel  half 
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width  h,  friction  velocity  uT,  wall  temperature  and  bulk  (integrated)  density  are  the  normalizing 
quantities  for  non-dimensionalization  with  a Reynolds  number  of  180.  Grid  refinement  studies 
(three  grids)  were  conducted  and  comparisons  were  made  with  results  from  incompressible  flow 
spectral  method  calculations  using  the  same  computational  box  size  and  the  finest  of  the  three 
grids.  Very  good  agreement  was  obtained.  In  addition,  using  the  same  uniform  body  force  term, 
computational  box  and  grid  size  as  in  Coleman  et  al.,  excellent  agreement  was  also  obtained 
with  the  Coleman  et  al.  spectral  compressible  Navier-Stokes  computation.  What  is  interesting 
is  that  this  simulation  did  not  require  filtering,  upwinding,  or  additional  numerical  dissipation 
for  shock  free  compressible  turbulence  computations.  Moreover,  this  high  order  method  can  be 
efficiently  extended  to  general  geometries  [28].  Using  CEN44-RK4-ENT  (f3  = 4)  together  with 
a so  called  Laplacian  viscous  term  formulation,  Sandham  & Yee  [19]  demonstrated  a robustness 
down  to  very  coarse  resolutions,  comparable  with  the  best  incompressible  turbulent  flow  solvers 
incorporating  de-aliasing  and  skew-symmetric  formulation  of  the  convection  terms.  For  the  same 
3-D  problem,  WEN05  is  more  than  six  times  as  expensive  yet  more  diffusive  than  the  present 
scheme  using  the  same  temporal  discretization.  Without  the  use  of  the  entropy  splitting  of  the 
inviscid  flux  derivatives  and  Laplacian  right  hand  side  formulation,  using  the  same  CFL  number, 
the  CEN44-RK4  solutions  diverge  for  all  three  grids  before  meaningful  turbulence  statistics  can 
be  obtained.  In  view  of  the  spurious  high-frequency-oscillation-producing  nonlinear  instability 
nature  of  central  schemes  and  the  past  experience  in  incompressible  turbulence  simulations,  in 
the  absence  of  entropy  splitting  and  Laplacian  formulation,  the  present  calculations  might  not 
even  be  possible  with  much  higher  resolution  grid. 

For  the  performance  of  ACM66-RK4,  WAV66-RK4  and  WEN05-RK4  on  a spatially  or  a 
time-developing  mixing  layer  problem  containing  shock  waves,  see  [31,  21], 

5.3.  MULTISCALE  COMPLEX  UNSTEADY  VISCOUS  COMPRESSIBLE  FLOWS  [22,  24] 

Extensive  grid  convergence  studies  using  WAV66-RK4  and  ACM66-RK4  for  two  complex  highly 
uusteady  viscous  compressible  flows  are  given  in  [22,  24].  The  first  flow  is  a 2-D  complex  viscous 
shock/shear/boundary- layer  interaction.  This  is  the  same  problem  and  flow  conditions  studied 
in  Daru  & Tenaud  [5].  The  second  flow  is  a supersonic  viscous  reacting  flow  concerning  fuel 
breakup.  More  accurate  solutions  were  obtained  with  WAVG6-RK4  and  ACM66-RK4  than  with 
WEN05-RK4,  which  is  nearly  three  times  as  expensive.  To  illustrate  the  performance  of  these 
nonlinear  filter  schemes,  the  first  model  is  considered.  The  ideal  gas  compressible  full  Navier- 
Stokes  equations  with  no  slip  BCs  at  the  adiabatic  walls  are  used.  The  fluid  is  at  rest  in  a 2-D 
box  0 < x,  y < 1.  A membrane  with  a shock  Mach  number  of  2.37  located  at  x = 1/2  separates 
two  different  states  of  the  gas.  The  dimensionless  initial  states  are 

PL  — 120,  pi  = 120/7;  PR  — 1-2,  Pr  = L2/7,  (5.1) 

where  pi,pr,  are  the  density  and  pressure  respectively,  to  the  left  of  x = 1/2,  and  pr,Pr  are 
the  same  quantities  to  the  right  of  x = 1/2.  7 = 1.4  and  the  Prandtl  number  is  0.73.  The  two 
Reynolds  numbers  considered  are  200  and  1000.  The  viscosity  is  assumed  to  be  constant  and 
independent  of  temperature,  so  Sutherland’s  law  is  not  used.  The  velocities  and  the  normal 
derivative  of  the  temperature  at  the  boundaries  are  set  equal  to  zero.  This  is  done  by  leaving  the 
value  of  the  density  obtained  by  the  one  sided  difference  scheme  at  the  boundary  unchanged, 
and  updating  the  energy  at  the  boundary  to  make  the  temperature  derivative  equal  to  zero. 

At  time  zero  the  membrane  is  removed  and  wave  interaction  occurs.  An  expansion  wave  and 
a shock  are  formed  initially.  Then,  a boundary  layer  is  formed  on  the  lower  boundary  behind  the 
right  going  waves.  After  reflection,  the  left  going  shock  wave  interacts  with  the  newly  formed 
boundary  layer,  causing  a number  of  vortices  and  lambda  shocks  near  the  boundary  layer.  Other 
kinds  of  layers  remain  after  the  shock  reflection  near  the  right  wall.  The  complexity  of  this  highly 
unsteady  shock/shear/boundary-layer  interactions  increases  as  the  Reynolds  number  increases. 
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For  illustration,  here  we  show  the  difficult  case  of  Reynolds  number  Re  = 1000.  The  com- 
putations stop  at  the  dimensionless  time  1 when  the  reflected  shock  wave  has  almost  reached 
the  middle  of  the  domain,  x = 1/2.  The  numerical  results  discussed  here  are  at  time  1 with 
uniform  Cartesian  grid  spacings  as  described  by  Daru  and  Tenaud.  Due  to  symmetry,  only  the 
lower  half  of  the  domain  is  used  in  the  computations;  symmetry  BCs  are  enforced  at  the  bound- 
ary y = 1/2.  Figure  5.5  show's  the  comparison  of  a second-order  MUSCL  using  a second-order 
Runge-Kutta  method  (MUSCL-RK2)  with  WAV66-RK4,  ACM66-RK4  and  WEN05-RK4  using 
a 1000  x 500  grid.  Comparing  with  the  converged  solution  of  WAV66-RK4  and  ACM66-RK4 
using  3000  x 1500  (see  bottom  of  figure)  and  4000  x 2000  grids  (see  [22]),  one  can  conclude 
that  VVAV66-RK4  exhibits  the  most  accurate  result  among  the  1000  x 500  grid  computations. 
We  note  that,  for  this  Reynolds  number,  the  unsteady  problem  is  extremely  stiff,  requiring  very 
small  time  steps  and  very  long-time  integrations  before  reaching  the  dimensionless  time  of  1. 


5.4.  AN  ADAPTIVE  NUMERICAL  DISSIPATION  MODEL  FOR  1-D  SHOCK-TURBULENCE 
INTERACTIONS 

In  classical  CFD  codes,  a second  order  accurate  base  method  is  used  together  with  two  con- 
stant strength  linear  numerical  dissipation  terms.  One  linear  fourth-order  dissipation  is  used 
everywhere  except  near  shocks/shears/steep-gradients  to  remove  nonlinear  instabilities.  It  does 
not  affect  the  second  order  accuracy  of  the  base  scheme.  The  second  dissipation  term  is  a 
second-order  linear  dissipation,  which  affects  the  order  of  accuracy,  but  is  only  switched  on  near 
discontinuities,  and/or  steep  unresolved  gradients  using  a gradient  sensor.  The  sensor  used  can- 
not distinguish  the  different  flow'  features  distinctly  and  is  not  accurate  enough  for  turbulent 
statistics  and  long-time  acoustic  computations,  unless  extreme  grid  refinement  is  employed. 

In  analogy  with  the  aforementioned  classical  methods,  a more  advanced  numerical  dissipation 
model  with  improved  flow  feature  extraction  sensors  for  high  order  central  schemes  is  proposed. 
Here,  we  consider  a dissipation  model  with  two  parts.  One  part  is  a nonlinear  filter  ([30])  and  the 
second  part  is  a high  order  linear  numerical  dissipation  term  modified  at  boundaries  to  become 
a scini-bounded  operator,  see  [20,  23].  The  wavelet  dissipation  control  sensor  developed  in  [21] 
is  used  as  the  flow  feature  detector.  The  sensor  is  computed  from  the  wavelet  estimate  of  the 
Lipschitz  exponent  a of  the  density  and  pressure  in  the  flow  field.  Below  we  present  a filter 
model  for  1-D  shock/ turbulence  interactions. 

A Filter  Model:  Using  a suitable  wavelet  basis  function,  the  final  result  of  the  wavelet  com- 
putation is  a quantity,  Sj,  which  is  near  zero  at  points  Xj  where  the  solution  is  smooth  and  near 
one  w'here  the  solution  is  discontinuous.  Sj  depends  on  the  Lipschitz  regularity  exponent  of  the 
solution.  We  define  the  filter  numerical  flux  of  the  numerical  dissipation  operator  as 

= mzx.(Sj,Sj-i)F-_l/2  + dj[l  - ma x(Sj, Sj-i)](h6 D+(D+D-)3Uj,  (5.2) 

where  F*  . is  the  flux  function  corresponding  to  the  dissipative  portion  of  a shock-capturing 
3 ‘■f* 

scheme  (e.g.,  second  order  accurate  TVD  scheme)  [30].  The  first  part  of  the  filter  stabilizes 
the  scheme  at  shock/shear  locations.  The  second  part,  is  an  eighth-order  filter  which  improves 
nonlinear  stability  away  from  shock/shear  locations.  Analogous  eighth-order  filters  can  be  used 
if  a sixth-order  compact  spatial  base  scheme  is  used  [7,  27].  We  switch  on  the  high  order  part  of 
the  filter  when  we  switch  off  the  nonlinear  filter.  The  physical  quantity  (e.g.,  local  Mach  number) 
can  be  used  to  determine  the  dj  parameter  of  this  high  order  dissipation  term. 

To  further  increase  stability  properties,  it  is  possible  to  use  the  sensor  to  switch  on  and  off 
the  entropy  splitting  and  adjust  the  value  of  the  entropy  splitting  parameter  according  to  flow 
type  and  region.  For  this  problem,  however,  we  believe  a constant  ft  = 1 away  from  the  shock 
waves  is  sufficient.  After  the  completion  of  a full  time  step  computation  using  the  sixth-order 
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a satisfying 


sup 

hf-0 


I f(x  + h)  -f{x)  I 
ha 


<C, 


(5.4) 


and  this  gives  information  about  the  regularity  of  the  function  / where  small  a means  poor 
regularity.  For  a Cl  wavelet  function  i/>  with  compact  support,  a can  be  estimated  from  the 
wavelet  coefficients,  defined  as 


=<  Sli’mj  > = J f(x)lpm,j(x)dx 


(5.5) 


where 

VV,  = 2mV<^)  (5-6) 

is  the  wavelet  function  i on  scale  m located  at  the  point  j in  space.  This  definition  gives  a so 
called  redundant  wavelet,  which  gives  (under  a few  technical  assumptions  on  rp)  a non-orthogonal 
basis  for  L 2.  It  is  possible  to  prove  that  the  coefficients  max,  \v>m  j \ in  a neighborhood  of  j0  decay 
as  2mo  as  the  scale  is  refined,  where  a is  the  Lipschitz  exponent  at  jo-  In  practical  computation, 
we  have  a smallest  scale,  determined  by  the  grid  size.  We  evaluate  wm,j  on  this  scale,  mg,  and 
a few  coarser  scales,  mg  + l,mo  + 2,  and  estimate  the  Lipschitz  exponent  at  the  point  jo  by  a 
least  square  fit  to  the  line  [21] 


max  log2  | Wm,7 1 = rnotj0  + c.  (5.7) 

j near  jo 

For  the  numerical  experiments,  the  wavelet  coefficient  wm,j  is  computed  numerically  by  a re- 
cursive procedure,  which  is  a second-order  B-spline  wavelet  or  a modification  of  Harten’s  multi- 
resolution scheme  [21].  The  sensor  we  use  in  the  computations  is 


Sj  = r(nj), 


(5.8) 


where 


r(aj) 


1 aj  < 0.5 

0 cij  > 0.5. 


(5.9) 


The  dissipative  model  (5.2)  is  used  to  solve  a simple,  yet  difficult,  1-D  compressible  inviscid 
shock-turbulence  interaction  problem  with  initial  data  consisting  of  a shock  propagating  into  an 
oscillatory  density.  The  initial  data  are  given  by 


(PL,  uL,  pL)  = (3.857143,  2.629369,  10.33333)  (5.10) 

to  the  left  of  a shock  located  at  x — —4,  and 

(pH,  UR,  pH)  = (1  + 0.2  *sin(5*x),  0,  1)  (5.11) 


to  the  right  of  the  shock  where  u is  the  velocity.  Fig.  6.6  show  the  comparison  between  a 
second-order  MUSCL-RK2  with  a sixth-order  central  scheme  and  the  aforementioned  numerical 
dissipation  model  using  R.K4  as  the  time  discretization  (WAV66-RK4-D8).  The  parameter  d = 
0.002  is  scaled  with  the  spectral  radius  of  the  Jacobian  of  the  ilux  function.  Note  that  the  eighth- 
order  dissipation  is  a filter,  and  is  different  from  the  CEN66-D8  used  in  Section  5.1  where  the 
dissipation  is  part  of  the  scheme.  Solution  using  a second-order  uniformly  non-oscillatory  (UNO) 
scheme  on  a 4000  uniform  grid  is  used  as  the  reference  solution  (solid  lines  on  the  first  three  sub- 
figures). The  bottom  of  the  right  figures  show  the  density  and  Lipschitz  exponent  distribution 
for  the  WAV66-RK4-D8  using  400  grid  points.  Comparing  our  result  with  the  most  accurate 
computation  found  in  the  literature  for  this  problem,  the  current  approach  is  highly  efficient 
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and  accurate  using  only  800  grid  points  without  grid  adaptation  or  very  high  order  shock- 
capturing scheme.  For  the  present  computation,  the  WAV66-RK4-D8  consumed  only  slightly 
more  CPU  than  the  second-order  scheme  MUSCL-RK2.  With  the  eighth-order  dissipation  filter 
turned  off  (i.e.,  only  the  nonlinear  filter  remains  - WAV66-RK4),  the  computation  is  not  very 
stable  unless  a finer  grid  and  smaller  time  step  is  used.  Turning  on  the  entropy  splitting  away 
from  the  shocks  helps  to  reduce  the  amount  of  the  eighth-order  dissipation  coefficient  [33]. 

6.  Concluding  Remarks 

A general  framework  for  the  design  of  an  adaptive  low  dissipative  high  order  scheme  is  pre- 
sented. The  approach  is  applicable  to  a wide  spectrum  of  flow  problems.  However  the  demand 
on  the  overall  numerical  approach  for  nonlinear  stability  and  accuracy  is  much  more  stringent 
for  long-time  integration  of  complex  multiscalc  shock/shcar/turbulcncc/acoustics  interactions 
and  numerical  combustion  problems.  Robust  classical  numerical  methods  for  less  complex  flow 
physics  are  not  suitable  or  practical  for  such  applications.  The  present  approach  is  designed 
expressly  to  address  such  flow  problems  and  computational  challenges.  The  incremental  stud- 
ies to  illustrate  the  performance  of  the  approach  are  summarized.  Extensive  testing  and  full 
implementation  of  the  approach  is  forthcoming.  The  results  shown  so  far  are  very  encouraging. 
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